
module module_CLM_HYDRO
    use domainMod        , only : ldomain
    use clmtype   , only : clm3
    use clm_varpar   , only : nlevgrnd
    use decompMod    , only : get_proc_bounds, get_proc_global
    use clm_varcon  , only : zsoi
    use subgridAveMod, only: l2g_1d
    use subgridAveMod, only: c2g_1d, c2g_2d


! NDHMS  module
     use module_mpp_land, only: global_nx, global_ny, decompose_data_real, &
                 write_io_real, my_id
    use module_hydro_stop, only: HYDRO_stop
    use module_HYDRO_drv, only: HYDRO_ini, HYDRO_exe

    implicit none
    integer begg, endg
    integer :: numg, numl, numc, nump
    INTEGER, PARAMETER :: double=8
    real(kind=double), pointer :: r2p(:,:) , r1p(:)

    integer ::  begl, endl, begc, endc, begp, endp
    real(kind=double), allocatable, dimension(:) :: clm_g1d
    real(kind=double), allocatable, dimension(:,:) :: clm_g2d

    real, allocatable, dimension(:,:) :: vg_test



CONTAINS

    subroutine clm_cpl_HYDRO()

       use module_rt_data, only:  rt_domain
       use module_CPL_LAND, only: CPL_LAND_INIT, cpl_outdate
       use module_mpp_land
       use config_base, only: nlst

        implicit none
        integer k, ix,jx, nn

        integer ::  did

        integer ntime, wrf_ix,wrf_jx
        real cpl_land_dt

        integer :: i,j
        integer clm_lev

!output flux and state variable

        did = 1

        call get_cpl_date(cpl_land_dt)

#ifdef HYDRO_D
        write(6,*) "cpl_land_dt = ",cpl_land_dt
#endif

        ntime = 1

        if(.not. RT_DOMAIN(did)%initialized) then




            call HYDRO_ini(ntime,did,1,1)

            if(nlst(did)%sys_cpl .ne. 4) then
               write(6,*) "FATAL ERROR: sys_cpl should be 4."
               call hydro_stop("In module_clm_HYDRO.F clm_cpl_HYDRO() - "// &
                               "sys_cpl should be 4.  Check hydro.namelist file")
            endif

            RT_DOMAIN(did)%initialized = .true.
            nlst(did)%dt = cpl_land_dt
            nlst(did)%startdate(1:19) = cpl_outdate(1:19)
            nlst(did)%olddate(1:19) = cpl_outdate(1:19)
        endif

        if(nlst(did)%rtFlag .eq. 0) return

        ix = rt_domain(did)%ix
        jx = rt_domain(did)%jx

        ! get the initial data from CLM
            call get_proc_global(numg, numl, numc, nump)

            call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp)
#ifdef HYDRO_D
            write(6,*) "begg, endg, ", begg, endg
            write(6,*) "begl, endl, ", begl, endl
            write(6,*) "begc, endc, ", begc, endc
            write(6,*) "begp, endp, ", begp, endp
#endif

!            call get_proc_bounds( begg=begg, endg=endg )
            nn = endg - begg + 1

            allocate(clm_g1d(endg-begg + 1))
            allocate(clm_g2d(endg-begg+1, nlevgrnd) )





            r1p => clm3%g%l%c%cwf%qflx_surf
            call c2g_1d(begc, endc, begl, endl, begg, endg, r1p(begc:endc), clm_g1d, &
                'urbanf', 'unity')
            call clm2ND2d(clm_g1d,nn,RT_DOMAIN(did)%INFXSRT,ix,jx)

#ifdef HYDRO_D
            write(6,*) "finish qflx_surf"

#endif

            r1p => clm3%g%l%c%cwf%qflx_drain
            call c2g_1d(begc, endc, begl, endl, begg, endg, r1p(begc:endc), clm_g1d, &
                'urbanf', 'unity')
            call clm2ND2d(clm_g1d,nn,RT_DOMAIN(did)%SOLDRAIN,ix,jx)


#ifdef HYDRO_D
            write(6,*) "finish qflx_drain "
#endif
            r2p =>clm3%g%l%c%ces%t_soisno
                call c2g_2d(begc, endc, begl, endl, begg, endg, nlevgrnd, r2p(begc:endc,1:nlevgrnd), clm_g2d, &
                'urbanh', 'unity')
            call clm2ND3d (clm_g2d,nn,nlevgrnd,zsoi(1:nlevgrnd), &
                      ix,jx, nlst(did)%nsoil,abs(rt_domain(did)%subsurface%properties%zsoil(1:nlst(did)%nsoil)),&
                      RT_DOMAIN(did)%STC)

#ifdef HYDRO_D
            write(6,*) "finish t_soisno   "
#endif

            r2p => clm3%g%l%c%cws%h2osoi_vol
                call c2g_2d(begc, endc, begl, endl, begg, endg, nlevgrnd, r2p(begc:endc,1:nlevgrnd), clm_g2d, &
                'urbanh', 'unity')
            call clm2ND3d (clm_g2d,nn,nlevgrnd,zsoi(1:nlevgrnd), &
                     ix,jx, nlst(did)%nsoil,abs(rt_domain(did)%subsurface%properties%zsoil(1:nlst(did)%nsoil)),&
                     RT_DOMAIN(did)%SMC)
#ifdef HYDRO_D
            write(6,*) "finish h2osoi_vol "

#endif

            r2p => clm3%g%l%c%cws%h2osoi_liq
            call c2g_2d(begc, endc, begl, endl, begg, endg, nlevgrnd, r2p(begc:endc,1:nlevgrnd), clm_g2d, &
                'urbanh', 'unity')
            call clm2ND3d (clm_g2d, nn,nlevgrnd,zsoi(1:nlevgrnd), &
                     ix,jx, nlst(did)%nsoil,abs(rt_domain(did)%subsurface%properties%zsoil(1:nlst(did)%nsoil)),&
                     RT_DOMAIN(did)%SH2OX)
            RT_DOMAIN(did)%SH2OX =  RT_DOMAIN(did)%SH2OX /1000.
 #ifdef HYDRO_D
            write(6,*) "finish h2osoi_liq "

            write(6,*) "before call HYDRO_exe"

#endif


            call HYDRO_exe(did)

#ifdef HYDRO_D
            write(6,*) "after call HYDRO_exe"

#endif
! add for update the CLM state variable.

           clm_lev = 10

            r2p =>clm3%g%l%c%ces%t_soisno
            call Toclm3d (clm_g2d,nn,clm_lev,zsoi(1:clm_lev), &
                     ix,jx, nlst(did)%nsoil,abs(rt_domain(did)%subsurface%properties%zsoil(1:nlst(did)%nsoil)),&
                     RT_DOMAIN(did)%STC)
            call g2c_2d(begc, endc, begl, endl, begg, endg,clm_lev,r2p(begc:endc,1:clm_lev), clm_g2d, &
                'urbanh', 'unity')




            r2p => clm3%g%l%c%cws%h2osoi_vol
            call Toclm3d (clm_g2d,nn,clm_lev,zsoi(1:clm_lev), &
                     ix,jx, nlst(did)%nsoil,abs(rt_domain(did)%subsurface%properties%zsoil(1:nlst(did)%nsoil)), &
                     RT_DOMAIN(did)%SMC)
            call g2c_2d_tmp(begc, endc, begl, endl, begg, endg,clm_lev,r2p(begc:endc,1:clm_lev), clm_g2d, &
                'urbanh', 'unity', 0.01, 10.)



            r2p => clm3%g%l%c%cws%h2osoi_liq
            call Toclm3d(clm_g2d,nn,clm_lev,zsoi(1:clm_lev), &
                     ix,jx, nlst(did)%nsoil,abs(rt_domain(did)%subsurface%properties%zsoil(1:nlst(did)%nsoil)),&
                     RT_DOMAIN(did)%SH2OX*1000)
            call g2c_2d_tmp(begc, endc, begl, endl, begg, endg,clm_lev,r2p(begc:endc,1:clm_lev), clm_g2d, &
                'urbanh', 'unity',0.01, 10.)




     ! 2d variable
     !      t_soisno
     ! 3 d variable
     deallocate(clm_g1d)
     deallocate(clm_g2d)
#ifdef HYDRO_D
            write(6,*) "end of drive ndhms"
#endif

     end subroutine clm_cpl_HYDRO

      subroutine Toclm3d (v1d_out,nn,kk1,z1_in,ix,jx,kk2,z2,v2_in)
        implicit none
          integer :: nn,ix,jx, kk2, kk1
          integer :: k
          real:: v2_in(ix,jx,kk2)
          real:: v1(ix,jx,kk1)
          real, dimension(kk2) :: z2
          REAL (KIND=double),dimension(nn,kk1):: v1d_out
          REAL (KIND=double),dimension(kk1):: z1_in
          REAL ,dimension(kk1):: z1

          z1 = z1_in

! v1 is the out from interp3D
             call Interp3D (z2,v2_in,kk2,z1,v1,ix,jx,kk1)



             do k = 1, kk1
                call Toclm2d(v1(:,:,k),ix,jx,v1d_out(:,k),nn)
             end do

      end subroutine Toclm3d

      subroutine Toclm2d(v1,ix,jx,v1d_out,nn)
         use module_mpp_land, only: my_id,IO_id,mpp_land_bcast_real
         implicit none
         integer ix,jx, nn
         real v1(ix,jx)
         REAL (KIND=double),dimension(nn) :: v1d_out
         real, dimension(global_nx*global_ny)::vg
         real, dimension(global_nx,global_ny)::vg2d


            call write_io_real(v1,vg2d)

            vg = 0.0
            if(my_id .eq. IO_id) then
                call TO1d(vg2d,vg,global_nx,global_ny)
            end if

#ifdef HYDRO_D
            write(6,*) "before scatter_1d_r"

#endif

            call scatter_1d_r(v1d_out,nn,vg,global_nx*global_ny)

#ifdef HYDRO_D
            write(6,*) "after scatter_1d_r"

#endif
      end subroutine Toclm2d

      subroutine TO1d(vg2d,v1d,nx,ny)
         implicit none
         integer nx,ny, n, i,j
         real vg2d(nx,ny)
         real v1d(nx*ny)
             do j = 1,ny
             do i = 1,nx
                n = (j-1)*nx + i
                v1d(n) = vg2d(i,j)
             enddo
             enddo
      end subroutine TO1d

      subroutine clm2ND3d (v1d_in,nn,kk1,z1_in,ix,jx,kk2,z2,vout)
        implicit none
          integer :: nn,kk1,ix,jx,kk2
          integer :: k
          real:: vout(ix,jx,kk2)
          real:: v1(ix,jx,kk1)
          real(kind=double),dimension(nn,kk1):: v1d_in
          real(kind=double), dimension(kk1) :: z1_in
          real, dimension(kk1) :: z1
          real, dimension(kk2) :: z2


          z1 = z1_in

          do k = 1, kk1
             call clm2ND2d(v1d_in(:,k),nn,v1(:,:,k),ix,jx)
          end do
             call Interp3D (z1(1:kk1),v1,kk1,z2(1:kk2),vout,ix,jx,kk2)

      end subroutine clm2ND3d

      subroutine clm2ND2d (v1d_in,nn,vout,ix,jx)
          use clmtype, only: grlnd
          use module_mpp_land, only: my_id,IO_id,mpp_land_bcast_real
          implicit none
          integer :: ix,jx,nn, n
          real vout(ix,jx)
          integer :: i,j,gni,gnj
!         ! real, allocatable, dimension(:)  :: g1d
          real(kind=double), dimension(global_nx*global_ny)  :: g1d
          real v2d_g(global_nx,global_ny)
          real(kind=double),dimension(nn) ::  v1d_in


          gni = global_nx
          gnj = global_ny
!     gather 1d global array

          if(endg .gt. gni*gnj) then
            write(6,*) "WARNING: endg > gni*gnj", endg, gni*gnj
            ! need to stop
          endif
          call gather_1d_real(v1d_in,nn,g1d,gni*gnj)
!  transfer 1d   to 2d

          if(my_id.eq.IO_id) then
             do j = 1,gnj
             do i = 1,gni
                n = (j-1)*gni + i
                v2d_g(i,j) = g1d(n)
             enddo
             enddo
          endif

!         if(my_id .eq. 0) then
!            call output_nc(v2d_g,gni,gnj, "test", "test.nc")
!         endif

! domain decomposition
          call decompose_data_real(v2d_g,vout)

      end subroutine clm2ND2d

      subroutine Interp3D (z1,v1,kk1,z,vout,ix,jx,kk)
!  input: z1,v1,kk1,z,ix,jx,kk
!  output: vout
!  interpolate based on soil layer: z1 and z
!  z :  soil layer of output variable.
!  z1: array of soil layers of input variable.
         implicit none
         integer:: i,j,k
         integer:: kk1, ix,jx,kk
         real :: z1(kk1), z(kk), v1(ix,jx,kk1),vout(ix,jx,kk)

!        write(6,*) "k, kk1 ", k, kk1
!        write(6,*) "z1(kk1), z(k) ", z1(kk1), z(k)

         do j = 1, jx
            do i = 1, ix
              do k = 1, kk
!                call interpLayer(abs(Z1),v1(i,j,1:kk1),kk1,abs( Z(k) ),vout(i,j,k))
                call interpLayer(Z1(1:kk1),v1(i,j,1:kk1),kk1,Z(k),vout(i,j,k))
              end do
            end do
         end do
      end subroutine Interp3D

      subroutine interpLayer(inZ,inV,inK,outZ,outV)
         implicit none
         integer:: k, k1, k2
         integer :: inK
         real:: inV(inK),inZ(inK)
         real:: outV, outZ, w1, w2

         if(outZ .le. inZ(1)) then
             if(inZ(2) .eq. inZ(1)) then
                write(6,*) "FATAL ERROR: inZ(2)=inZ(1) ", inZ(2),inZ(1)
               !stop 99
               stop("In module_clm_HYDRO.F interpLayer() - inZ(2)=inZ(1)")
             end if
             w1 = (inZ(2)-outZ)/(inZ(2)-inZ(1))
             w2 = (inZ(1)-outZ)/(inZ(2)-inZ(1))
             outV = inV(1)*w1-inV(2)*w2
             return
         elseif(outZ .ge. inZ(inK)) then
             if(inZ(inK-1) .eq. inZ(inK)) then
                write(6,*) "FATAL ERROR: inZ(inK-1)=inZ(inK) ", inZ(inK-1),inZ(inK)
                !stop 99
                stop("FATAL ERROR: In module_clm_HYDRO.F interpLayer() - inZ(inK-1)=inZ(inK)")
             end if
             w1 = (outZ-inZ(inK-1))/(inZ(inK)-inZ(inK-1))
             w2 = (outZ-inZ(inK))  /(inZ(inK)-inZ(inK-1))
             outV = inV(inK)*w1 -inV(inK-1)* w2
             return
         else
            do k = 2, inK
             if((inZ(k) .ge. outZ).and.(inZ(k-1) .le. outZ) ) then
                k1  = k-1
                k2 = k
             if(inZ(k1) .eq. inZ(k2)) then
                write(6,*) "FATAL ERROR: inZ(k1)=inZ(k2) ", inZ(k1),inZ(k2)
                !stop 99
                stop("FATAL ERROR: In module_clm_HYDRO.F interpLayer()- inZ(k1)=inZ(k2)")
             end if
                w1 = (outZ-inZ(k1))/(inZ(k2)-inZ(k1))
                w2 = (inZ(k2)-outZ)/(inZ(k2)-inZ(k1))
                outV = inV(k2)*w1 + inV(k1)*w2
                return
             end if
            end do
         endif
      end subroutine interpLayer

    subroutine get_cpl_date(cpl_land_dt)
      use clm_time_manager, only : get_step_size, get_curr_date
      use module_CPL_LAND, only: cpl_outdate
      implicit none
      integer :: yr, mo, da, hr, mn, ss , sec
      real :: cpl_land_dt
      integer cmdate

      cpl_land_dt  = 1.0* get_step_size()

      call get_curr_date(yr, mo, da, sec)


      hr = sec/3600

      sec = mod(sec,3600)
      mn = sec/60

      ss = mod(sec,60)

      if(yr .lt. 1000) yr = 2000+yr
      write(cpl_outdate(1:4),'(I4)') yr

      if(mo .lt. 10 ) then
         write(cpl_outdate(5:7),'("-0",I1)') mo
      else
         write(cpl_outdate(5:7),'("-",I2)') mo
      endif

      if(da .lt. 10 ) then
         write(cpl_outdate(8:10),'("-0",I1)') da
      else
         write(cpl_outdate(8:10),'("-",I2)') da
      endif

      if(hr .lt. 10 ) then
         write(cpl_outdate(11:13),'("_0",I1)') hr
      else
         write(cpl_outdate(11:13),'("_",I2)') hr
      endif

      if(mn .lt. 10 ) then
         write(cpl_outdate(14:16),'(":0",I1)') mn
      else
         write(cpl_outdate(14:16),'(":",I2)') mn
      endif

      if(ss .lt. 10 ) then
         write(cpl_outdate(17:19),'(":0",I1)') ss
      else
         write(cpl_outdate(17:19),'(":",I2)') ss
      endif



!      write(cpl_outdate,'(I4,"-",I2,"-",I2,"_",I2,":",I2,":",I2)' )     &
!            yr, mo, da, hr, mn, ss

#ifdef HYDRO_D
      write(6,*) "cpl_outdate = ",cpl_outdate

#endif
    end subroutine get_cpl_date

!    transfer the grid column to gridcell.
  subroutine g2c_1d(lbc, ubc, lbl, ubl, lbg, ubg, carr, garr, &
       c2l_scale_type, l2g_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from columns to gridcells.
! Averaging is only done for points that are not equal to "spval".
!
! !ARGUMENTS:
     use clm_varcon, only : spval, isturb,  icol_roof, icol_sunwall, icol_shadewall, &
                         icol_road_perv, icol_road_imperv
     use clm_varctl, only : iulog


    implicit none
    integer , intent(in)  :: lbc, ubc              ! beginning and ending column indices
    integer , intent(in)  :: lbl, ubl              ! beginning and ending landunit indices
    integer , intent(in)  :: lbg, ubg              ! beginning and ending landunit indices
    real(kind=double), intent(out)  :: carr(lbc:ubc)         ! output column array
    real(kind=double), intent(in) :: garr(lbg:ubg)         ! input  gridcell array
    character(len=*), intent(in) :: c2l_scale_type ! scale factor type for averaging
    character(len=*), intent(in) :: l2g_scale_type ! scale factor type for averaging
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: ci,c,l,g,index         ! indices
    integer  :: max_col_per_gcell      ! max columns per gridcell; on the fly
    logical  :: found                  ! temporary for error check
    real(kind=double) :: scale_c2l(lbc:ubc)     ! scale factor
    real(kind=double) :: scale_l2g(lbl:ubl)     ! scale factor
    real(kind=double) :: sumwt(lbg:ubg)         ! sum of weights
    real(kind=double), pointer :: wtgcell(:)    ! weight of columns relative to gridcells
    integer , pointer :: clandunit(:)  ! landunit of corresponding column
    integer , pointer :: cgridcell(:)  ! gridcell of corresponding column
    integer , pointer :: ncolumns(:)   ! number of columns in gridcell
    integer , pointer :: coli(:)       ! initial column index in gridcell
    integer , pointer :: ctype(:)      ! column type
    integer , pointer :: ltype(:)      ! landunit type
    real(kind=double), pointer :: canyon_hwr(:) ! urban canyon height to width ratio
    integer :: gcount(lbg:ubg)
!------------------------------------------------------------------------


    ctype      => clm3%g%l%c%itype
    ltype      => clm3%g%l%itype
    canyon_hwr => clm3%g%l%canyon_hwr
    wtgcell    => clm3%g%l%c%wtgcell
    clandunit  => clm3%g%l%c%landunit
    cgridcell  => clm3%g%l%c%gridcell
    ncolumns   => clm3%g%ncolumns
    coli       => clm3%g%coli

    if (l2g_scale_type == 'unity') then
       do l = lbl,ubl
          scale_l2g(l) = 1.0
       end do
    end if

    if (c2l_scale_type == 'unity') then
       do c = lbc,ubc
          scale_c2l(c) = 1.0
       end do
    else if (c2l_scale_type == 'urbanf') then
       do c = lbc,ubc
          l = clandunit(c)
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0
             end if
          else
             scale_c2l(c) = 1.0
          end if
       end do
    else if (c2l_scale_type == 'urbans') then
       do c = lbc,ubc
          l = clandunit(c)
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0 / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0
             end if
          else
             scale_c2l(c) = 1.0
          end if
       end do
    else if (c2l_scale_type == 'urbanh') then
       do c = lbc,ubc
          l = clandunit(c)
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = spval
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = spval
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = spval
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = spval
             end if
          else
             scale_c2l(c) = 1.0
          end if
       end do
    end if

    sumwt(:) = 0.
    gcount(:) = 0
    do c = lbc,ubc
          if (carr(c) /= spval .and. scale_c2l(c) /= spval) then
             l = clandunit(c)
             g = cgridcell(c)
             if(wtgcell(c) .ne. 0) then
                carr(c) = garr(g)
             endif
          end if
    end do
    end subroutine g2c_1d

    subroutine gather_1d_real(l1d,lnn,g1d,gnn)
       use spmdGathScatMod
       use clmtype   , only : grlnd, nameg

       implicit none
       integer lnn, gnn
       real(kind=double),dimension(lnn) :: l1d
       real(kind=double),dimension(gnn) :: g1d
       real(kind=double), pointer :: lp_yw(:), gp_yw(:)

       allocate(lp_yw(lnn) )
       allocate(gp_yw(gnn) )
       lp_yw = l1d


       call gather_data_to_master(lp_yw,gp_yw,grlnd)

       g1d = gp_yw
       deallocate(lp_yw)
       deallocate(gp_yw)
    end subroutine gather_1d_real

    subroutine scatter_1d_r(l1d,lnn,g1d,gnn)
       use spmdGathScatMod, only: scatter_1darray_real
       use clmtype   , only : grlnd, nameg

       implicit none
       integer lnn, gnn
       real(kind=double),dimension(lnn) :: l1d
       real,dimension(gnn) :: g1d
       real(kind=double), pointer :: lp_yw(:), gp_yw(:)




       allocate(lp_yw(begg:endg) )
       allocate(gp_yw(gnn) )


       gp_yw(:) = g1d(:)


!      call scatter_data_from_master(lp, gp, grlnd)

       call scatter_1darray_real(lp_yw, gp_yw, grlnd)


       l1d = lp_yw(begg:endg)

       deallocate(lp_yw)
       deallocate(gp_yw)
    end subroutine scatter_1d_r

      subroutine output_nc(array,idim,jdim, var_name, file_name)
      implicit none
#include <netcdf.inc>
          integer idim,jdim
          real array(idim,jdim)
          integer dim(2)
          character(len=*) file_name,var_name
          integer   iret,ncid,varid,idim_id,jdim_id
          integer i,j
          iret = nf_create(trim(file_name), 0, ncid)
          iret = nf_def_dim(ncid, "idim", idim,idim_id)
          iret = nf_def_dim(ncid, "jdim", jdim,jdim_id)
          dim(1)=idim_id
          dim(2)=jdim_id
          iret = nf_put_att_real(ncid, NF_GLOBAL,  &
                "missing_value", NF_FLOAT, 1, -1.E33)

          iret = nf_def_var(ncid,var_name,NF_FLOAT,2,dim,varid)
          iret = nf_enddef(ncid)

!   output
          iret = nf_inq_varid(ncid,var_name,varid)
          iret = nf_put_var_real(ncid,varid,array)
          iret=nf_close(ncid)
      end subroutine output_nc

  subroutine g2c_2d(lbc, ubc, lbl, ubl, lbg, ubg, num2d, carr, garr, &
       c2l_scale_type, l2g_scale_type)
!
! !DESCRIPTION:
! Perfrom subgrid-average from columns to gridcells.
! Averaging is only done for points that are not equal to "spval".
!
     use clm_varcon, only : spval, isturb,  icol_roof, icol_sunwall, icol_shadewall, &
                         icol_road_perv, icol_road_imperv
     use clm_varctl, only : iulog
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbc, ubc              ! beginning and ending column indices
    integer , intent(in)  :: lbl, ubl              ! beginning and ending landunit indices
    integer , intent(in)  :: lbg, ubg              ! beginning and ending gridcell indices
    integer , intent(in)  :: num2d                 ! size of second dimension
    real(kind=double), intent(inout)  :: carr(lbc:ubc,num2d)   ! input column array
    real(kind=double), intent(in) :: garr(lbg:ubg,num2d)   ! output gridcell array
    character(len=*), intent(in) :: c2l_scale_type ! scale factor type for averaging
    character(len=*), intent(in) :: l2g_scale_type ! scale factor type for averaging

    real(kind=double) :: yw_r
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: j,ci,c,g,l,index       ! indices
    integer  :: max_col_per_gcell      ! max columns per gridcell; on the fly
    logical  :: found                  ! temporary for error check
    real(kind=double) :: scale_c2l(lbc:ubc)     ! scale factor
    real(kind=double) :: scale_l2g(lbl:ubl)     ! scale factor
    real(kind=double) :: sumwt(lbg:ubg)         ! sum of weights
    real(kind=double), pointer :: wtgcell(:)    ! weight of columns relative to gridcells
    integer , pointer :: clandunit(:)  ! landunit of corresponding column
    integer , pointer :: cgridcell(:)  ! gridcell of corresponding column
    integer , pointer :: ncolumns(:)   ! number of columns in gridcell
    integer , pointer :: coli(:)       ! initial column index in gridcell
    integer , pointer :: ctype(:)      ! column type
    integer , pointer :: ltype(:)      ! landunit type
    real(kind=double), pointer :: canyon_hwr(:) ! urban canyon height to width ratio
    real :: w_yw
!------------------------------------------------------------------------

    ctype      => clm3%g%l%c%itype
    ltype      => clm3%g%l%itype
    canyon_hwr => clm3%g%l%canyon_hwr
    wtgcell    => clm3%g%l%c%wtgcell
    clandunit  => clm3%g%l%c%landunit
    cgridcell  => clm3%g%l%c%gridcell
    ncolumns   => clm3%g%ncolumns
    coli       => clm3%g%coli

    if (l2g_scale_type == 'unity') then
       do l = lbl,ubl
          scale_l2g(l) = 1.0
       end do
    else
       write(iulog,*)'WARNING: In module_clm_HYDRO.F c2g_2d() - '// &
                     'scale type ',l2g_scale_type,' not supported'
    end if
    if (c2l_scale_type == 'unity') then
       do c = lbc,ubc
          scale_c2l(c) = 1.0
       end do
    else if (c2l_scale_type == 'urbanf') then
       do c = lbc,ubc
          l = clandunit(c)
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0
             end if
          else
             scale_c2l(c) = 1.0
          end if
       end do
    else if (c2l_scale_type == 'urbans') then
       do c = lbc,ubc
          l = clandunit(c)
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0 / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0
             end if
          else
             scale_c2l(c) = 1.0
          end if
       end do
    else if (c2l_scale_type == 'urbanh') then
       do c = lbc,ubc
          l = clandunit(c)
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = spval
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = spval
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = spval
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = spval
             end if
          else
             scale_c2l(c) = 1.0
          end if
       end do
    else
      write(iulog,*) 'WARNING: In module_clm_HYDRO.F c2g_2d()- '// &
                     'scale type ',c2l_scale_type,' not supported'
    end if

    do j = 1,num2d
       sumwt(:) = 0.
       do c = lbc,ubc
          if (wtgcell(c) /= 0) then
             if (carr(c,j) /= spval .and. scale_c2l(c) /= spval) then
                l = clandunit(c)
                g = cgridcell(c)
                sumwt(g) = sumwt(g) + wtgcell(c)
             end if
          end if
       end do
       do c = lbc,ubc
          if (wtgcell(c) /= 0) then
             if (carr(c,j) /= spval .and. scale_c2l(c) /= spval) then
                l = clandunit(c)
                g = cgridcell(c)
                   !garr(g,j) = garr(g,j) + carr(c,j) * scale_c2l(c) * scale_l2g(l) * wtgcell(c)
                   !w_yw = scale_c2l(c) * scale_l2g(l) * wtgcell(c)
                   w_yw = scale_c2l(c) * scale_l2g(l)
                   if(w_yw .ne. 0) then
                      ! carr(c,j) =  garr(g,j) / w_yw *sumwt(g)
                      yw_r = garr(g,j) / w_yw
                      if(abs(yw_r - carr(c,j)) .lt. 400 .and. yw_r .ne. 0 ) then
                         carr(c,j) =  garr(g,j) / w_yw
                      endif
                   endif
            end if
          end if
       end do
    end do

  end subroutine g2c_2d

  subroutine g2c_2d_tmp(lbc, ubc, lbl, ubl, lbg, ubg, num2d, carr, garr, &
       c2l_scale_type, l2g_scale_type, minv, maxv)
!
! !DESCRIPTION:
! Perfrom subgrid-average from columns to gridcells.
! Averaging is only done for points that are not equal to "spval".
!
     use clm_varcon, only : spval, isturb,  icol_roof, icol_sunwall, icol_shadewall, &
                         icol_road_perv, icol_road_imperv
     use clm_varctl, only : iulog
! !ARGUMENTS:
    implicit none
    integer , intent(in)  :: lbc, ubc              ! beginning and ending column indices
    integer , intent(in)  :: lbl, ubl              ! beginning and ending landunit indices
    integer , intent(in)  :: lbg, ubg              ! beginning and ending gridcell indices
    integer , intent(in)  :: num2d                 ! size of second dimension
    real(kind=double), intent(inout)  :: carr(lbc:ubc,num2d)   ! input column array
    real(kind=double), intent(in) :: garr(lbg:ubg,num2d)   ! output gridcell array
    character(len=*), intent(in) :: c2l_scale_type ! scale factor type for averaging
    character(len=*), intent(in) :: l2g_scale_type ! scale factor type for averaging

    real(kind=double) :: yw_r
    real :: minv, maxv
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 12/03
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: j,ci,c,g,l,index       ! indices
    integer  :: max_col_per_gcell      ! max columns per gridcell; on the fly
    logical  :: found                  ! temporary for error check
    real(kind=double) :: scale_c2l(lbc:ubc)     ! scale factor
    real(kind=double) :: scale_l2g(lbl:ubl)     ! scale factor
    real(kind=double) :: sumwt(lbg:ubg)         ! sum of weights
    real(kind=double), pointer :: wtgcell(:)    ! weight of columns relative to gridcells
    integer , pointer :: clandunit(:)  ! landunit of corresponding column
    integer , pointer :: cgridcell(:)  ! gridcell of corresponding column
    integer , pointer :: ncolumns(:)   ! number of columns in gridcell
    integer , pointer :: coli(:)       ! initial column index in gridcell
    integer , pointer :: ctype(:)      ! column type
    integer , pointer :: ltype(:)      ! landunit type
    real(kind=double), pointer :: canyon_hwr(:) ! urban canyon height to width ratio
    real :: w_yw
!------------------------------------------------------------------------

    ctype      => clm3%g%l%c%itype
    ltype      => clm3%g%l%itype
    canyon_hwr => clm3%g%l%canyon_hwr
    wtgcell    => clm3%g%l%c%wtgcell
    clandunit  => clm3%g%l%c%landunit
    cgridcell  => clm3%g%l%c%gridcell
    ncolumns   => clm3%g%ncolumns
    coli       => clm3%g%coli

    if (l2g_scale_type == 'unity') then
       do l = lbl,ubl
          scale_l2g(l) = 1.0
       end do
    else
      write(iulog,*) 'WARNING: In module_clm_HYDRO.F g2c_2d_tmp - '// &
                     'scale type ',l2g_scale_type,' not supported'
    end if
    if (c2l_scale_type == 'unity') then
       do c = lbc,ubc
          scale_c2l(c) = 1.0
       end do
    else if (c2l_scale_type == 'urbanf') then
       do c = lbc,ubc
          l = clandunit(c)
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = 3.0 * canyon_hwr(l)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0
             end if
          else
             scale_c2l(c) = 1.0
          end if
       end do
    else if (c2l_scale_type == 'urbans') then
       do c = lbc,ubc
          l = clandunit(c)
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = (3.0 * canyon_hwr(l)) / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = 3.0 / (2.*canyon_hwr(l) + 1.)
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = 1.0
             end if
          else
             scale_c2l(c) = 1.0
          end if
       end do
    else if (c2l_scale_type == 'urbanh') then
       do c = lbc,ubc
          l = clandunit(c)
          if (ltype(l) == isturb) then
             if (ctype(c) == icol_sunwall) then
                scale_c2l(c) = spval
             else if (ctype(c) == icol_shadewall) then
                scale_c2l(c) = spval
             else if (ctype(c) == icol_road_perv .or. ctype(c) == icol_road_imperv) then
                scale_c2l(c) = spval
             else if (ctype(c) == icol_roof) then
                scale_c2l(c) = spval
             end if
          else
             scale_c2l(c) = 1.0
          end if
       end do
    else
      write(iulog,*) 'WARNING: In module_clm_HYDRO.F g2c_2d_tmp() - '// &
                      'scale type ',c2l_scale_type,' not supported'
    end if

    do j = 1,num2d
       sumwt(:) = 0.
       do c = lbc,ubc
          if (wtgcell(c) /= 0) then
             if (carr(c,j) /= spval .and. scale_c2l(c) /= spval) then
                l = clandunit(c)
                g = cgridcell(c)
                sumwt(g) = sumwt(g) + wtgcell(c)
             end if
          end if
       end do
       do c = lbc,ubc
          if (wtgcell(c) /= 0) then
             if (carr(c,j) /= spval .and. scale_c2l(c) /= spval) then
                    l = clandunit(c)
                    g = cgridcell(c)
                    !garr(g,j) = garr(g,j) + carr(c,j) * scale_c2l(c) * scale_l2g(l) * wtgcell(c)
                    !w_yw = scale_c2l(c) * scale_l2g(l) * wtgcell(c)
                    w_yw = scale_c2l(c) * scale_l2g(l)
                    if(w_yw .ne. 0) then
                    if(abs(garr(g,j)) .gt. minv .and. abs(garr(g,j)) .lt. maxv) then
                       ! carr(c,j) =  garr(g,j) / w_yw *sumwt(g)
                             yw_r =  garr(g,j) / w_yw
                             if(yw_r .gt. 0) then
                                carr(c,j) =  garr(g,j) / w_yw
                             endif
                    endif
                    endif
             end if
          end if
       end do
    end do

  end subroutine g2c_2d_tmp
end module module_clm_HYDRO
